sortLvlsByVar.fnc <- function(oldFactor, sortingVariable, ascending = TRUE) {
  
  require("dplyr")
  require("magrittr")
  
  # Combine into data frame
  df <- data.frame(oldFactor, sortingVariable)
  
  ###
  ### If you want to sort the levels by, say, the median, sd etc. instead of the mean,
  ### just change 'mean(sortingVariable)' below to, say, 'median(sortingVariable)'.
  ###
  
  # Compute average of sortingVariable and arrange (ascending)
  if (ascending == TRUE) {
    df_av <- df %>% group_by(oldFactor) %>% summarise(meanSortingVariable = mean(sortingVariable)) %>% 
      arrange(meanSortingVariable)
  }
  
  # Compute average of sortingVariable and arrange (descending)
  if (ascending == FALSE) {
    df_av <- df %>% group_by(oldFactor) %>% summarise(meanSortingVariable = mean(sortingVariable)) %>% 
      arrange(desc(meanSortingVariable))
  }
  
  # Return factor with new level order
  newFactor <- factor(oldFactor, levels = df_av$oldFactor)
  return(newFactor)
}

# Sort factor levels by their frequency of occurrence
sortLvlsByN.fnc <- function(oldFactor, ascending = TRUE) {
  
  require("magrittr")
  
  # Return factor with new level order
  newFactor <- factor(oldFactor, levels = table(oldFactor)  %>% sort(., decreasing = !ascending)  %>% names())
  return(newFactor)
}

# Sort factor levels arbitrarily
sortLvls.fnc <- function(oldFactor, levelOrder) {
  if(!is.factor(oldFactor)) stop("The variable you want to reorder isn't a factor.")
  
  if(!is.numeric(levelOrder)) stop("'order' should be a numeric vector.")
  
  if(max(levelOrder) > length(levels(oldFactor))) stop("The largest number in 'order' can't be larger than the number of levels in the factor.")
  
  if(length(levelOrder) > length(levels(oldFactor))) stop("You can't have more elements in 'order' than there are levels in the factor.")
  
  if(length(levelOrder) == length(levels(oldFactor))) {
    reorderedFactor <- factor(oldFactor, levels = levels(oldFactor)[levelOrder])
  }
  
  if(length(levelOrder) < length(levels(oldFactor))) {
    levelOrderAll <- c(levelOrder, (1:length(levels(oldFactor)))[-levelOrder])
    reorderedFactor <- factor(oldFactor, levels = levels(oldFactor)[levelOrderAll])
  }
  
  return(reorderedFactor)
}

#############################################################################
require(tidyr)
require(dplyr)

out_path = "D:/Research/global_refor_mac/model_runs/"

model_list = readRDS(paste0(out_path, "glob_model_run_stand.rds"))

# function to create datasets by continent and under different price regimes
calc_data_model <- function(model_list,
                            year_vec,
                            dollar_prices){
  
  
  for_skm_list = list()
  def_skm_list = list()
  def_emiss_list = list()
  for_refor_skm_list = list()
  gain_skm_list = list()
  gain_emiss_list = list()
  
  for(dl in 1:length(dollar_prices)){
    
    # calculate forest square kilometers 
    for_skm_df = model_list[[dl]][[1]] %>%
      gather(year, fc_skm, `2020`:`2050`) %>%
      group_by(continent, year) %>%
      summarize(fc_skm = sum(fc_skm, na.rm = TRUE))
    for_skm_df$dollar = dollar_prices[dl]
    for_skm_list[[dl]] = for_skm_df
    
    # calculate deforestation square kilometers
    def_skm_df = model_list[[dl]][[2]] %>%
      gather(year, def_skm, `2020`:`2050`) %>%
      group_by(continent, year) %>%
      summarize(def_skm = sum(def_skm, na.rm = TRUE))
    def_skm_df$dollar = dollar_prices[dl]
    def_skm_list[[dl]]  = def_skm_df
    
    # calculate deforestation square kilometers
    def_emiss_df = model_list[[dl]][[3]] %>%
      gather(year, def_emiss, `2020`:`2050`) %>%
      group_by(continent, year) %>%
      summarize(def_emiss = sum(def_emiss, na.rm = TRUE))
    def_emiss_df$dollar = dollar_prices[dl]
    def_emiss_list[[dl]] = def_emiss_df
    
    # calculate forest square kilometers 
    for_skm_refor_df = model_list[[dl]][[4]] %>%
      gather(year, fc_skm, `2020`:`2050`) %>%
      group_by(continent, year) %>%
      summarize(fc_skm = sum(fc_skm, na.rm = TRUE))
    for_skm_refor_df$dollar = dollar_prices[dl]
    for_refor_skm_list[[dl]] = for_skm_refor_df
    
    # calculate deforestation square kilometers
    gain_skm_df = model_list[[dl]][[5]] %>%
      gather(year, gain_skm, `2020`:`2050`) %>%
      group_by(continent, year) %>%
      summarize(gain_skm = sum(gain_skm, na.rm = TRUE))
    gain_skm_df$dollar = dollar_prices[dl]
    gain_skm_list[[dl]] = gain_skm_df
    
    # calculate deforestation square kilometers
    gain_emiss_df = model_list[[dl]][[6]] %>%
      gather(year, gain_emiss, `2020`:`2050`) %>%
      group_by(continent, year) %>%
      summarize(gain_emiss = sum(gain_emiss, na.rm = TRUE))
    gain_emiss_df$dollar = dollar_prices[dl]
    gain_emiss_list[[dl]] = gain_emiss_df
    
  }
  
  for_skm_list = bind_rows(for_skm_list)
  def_skm_list = bind_rows(def_skm_list)
  def_emiss_list = bind_rows(def_emiss_list)
  for_refor_skm_list = bind_rows(for_refor_skm_list)
  gain_skm_list = bind_rows(gain_skm_list)
  gain_emiss_list = bind_rows(gain_emiss_list)
  
  return(list(for_skm_list, def_skm_list,
              def_emiss_list, for_refor_skm_list, gain_skm_list,
              gain_emiss_list))
  
}

# function to create datasets by continent and under different price regimes
calc_data_model_cntry <- function(model_list,
                            year_vec,
                            dollar_prices){

  
  for_skm_list = list()
  def_skm_list = list()
  def_emiss_list = list()
  for_refor_skm_list = list()
  gain_skm_list = list()
  gain_emiss_list = list()
  
  for(dl in 1:length(dollar_prices)){
    
    # calculate forest square kilometers 
    for_skm_df = model_list[[dl]][[1]] %>%
      gather(year, fc_skm, `2020`:`2050`) %>%
      group_by(country, year) %>%
      summarize(fc_skm = sum(fc_skm, na.rm = TRUE))
    for_skm_df$dollar = dollar_prices[dl]
    for_skm_list[[dl]] = for_skm_df
    
    # calculate deforestation square kilometers
    def_skm_df = model_list[[dl]][[2]] %>%
      gather(year, def_skm, `2020`:`2050`) %>%
      group_by(country, year) %>%
      summarize(def_skm = sum(def_skm, na.rm = TRUE))
    def_skm_df$dollar = dollar_prices[dl]
    def_skm_list[[dl]]  = def_skm_df
    
    # calculate deforestation square kilometers
    def_emiss_df = model_list[[dl]][[3]] %>%
      gather(year, def_emiss, `2020`:`2050`) %>%
      group_by(country, year) %>%
      summarize(def_emiss = sum(def_emiss, na.rm = TRUE))
    def_emiss_df$dollar = dollar_prices[dl]
    def_emiss_list[[dl]] = def_emiss_df
    
    # calculate forest square kilometers 
    for_skm_refor_df = model_list[[dl]][[4]] %>%
      gather(year, fc_skm, `2020`:`2050`) %>%
      group_by(country, year) %>%
      summarize(fc_skm = sum(fc_skm, na.rm = TRUE))
    for_skm_refor_df$dollar = dollar_prices[dl]
    for_refor_skm_list[[dl]] = for_skm_refor_df
    
    # calculate deforestation square kilometers
    gain_skm_df = model_list[[dl]][[5]] %>%
      gather(year, gain_skm, `2020`:`2050`) %>%
      group_by(country, year) %>%
      summarize(gain_skm = sum(gain_skm, na.rm = TRUE))
    gain_skm_df$dollar = dollar_prices[dl]
    gain_skm_list[[dl]] = gain_skm_df
    
    # calculate deforestation square kilometers
    gain_emiss_df = model_list[[dl]][[6]] %>%
      gather(year, gain_emiss, `2020`:`2050`) %>%
      group_by(country, year) %>%
      summarize(gain_emiss = sum(gain_emiss, na.rm = TRUE))
    gain_emiss_df$dollar = dollar_prices[dl]
    gain_emiss_list[[dl]] = gain_emiss_df
    
  }
  
  for_skm_list = bind_rows(for_skm_list)
  def_skm_list = bind_rows(def_skm_list)
  def_emiss_list = bind_rows(def_emiss_list)
  for_refor_skm_list = bind_rows(for_refor_skm_list)
  gain_skm_list = bind_rows(gain_skm_list)
  gain_emiss_list = bind_rows(gain_emiss_list)
  
  return(list(for_skm_list, def_skm_list,
              def_emiss_list, for_refor_skm_list, gain_skm_list,
              gain_emiss_list))
  
}

# read in the csv files
hist_data = read.csv("D:/Research/global_refor_mac/historical_data/historical_data.csv",
                     stringsAsFactors = FALSE)
hist_data = hist_data[which(!is.na(hist_data$cont_var)), ]

### create standard output
standard_output <- calc_data_model(model_list = model_list,
                                   year_vec = seq(2020, 2050, 10),
                                   dollar_prices = c(0,5,seq(10,100,10)))

standard_output_cntry <- calc_data_model_cntry(model_list = model_list,
                                   year_vec = seq(2020, 2050, 10),
                                   dollar_prices = c(0,5,seq(10,100,10)))

##### 1. Create predicted business as usual deforestation
################################################################

graphics_folder_loc = "D:/Research/global_refor_mac/graphics/"

require(ggplot2)
require(scales)
require(RColorBrewer)

# match names 
match_df = data.frame("continent" = c("Africa",
                                      "Asia",
                                      "Latin America",
                                      "Pan-tropical"),
                      "name" = c("Africa", "Asia",
                                 "LatinAmerica", "global"),
                      stringsAsFactors = FALSE)

pred_def = standard_output[[2]] %>% filter(dollar == 0)
global_def = pred_def %>% group_by(year) %>% summarize(def_skm = sum(def_skm,
                                                                     na.rm = TRUE)) %>%
  mutate(continent = "global",
         dollar = 0)

hist_def = hist_data %>% group_by(cont_var) %>% summarize(def_skm = sum(def_skm,
                                                                        na.rm = TRUE))  %>%
  mutate(year = 2010,
         dollar = 0)
global_hist = hist_def %>% summarize(def_skm = sum(def_skm),
                                     year = mean(year),
                                     dollar = mean(dollar)) %>%
  mutate(cont_var = "global")

# create a historical data.frame 
hist_def = bind_rows(hist_def, global_hist)
pred_def = bind_rows(pred_def, global_def)
pred_def$year = as.numeric(pred_def$year)
colnames(hist_def)[1] = "continent"

# match the names 
hist_def$clean_name = paste0(match_df$continent[match(hist_def$continent, match_df$name)],
                             " - historical")
pred_def$clean_name = paste0(match_df$continent[match(pred_def$continent, match_df$name)],
                             " - projected")

# graph data.frame
graph_df = bind_rows(hist_def, pred_def)
graph_df$continent_name = match_df$continent[match(graph_df$continent, match_df$name)]
graph_df$continent_name_fac = factor(graph_df$continent_name)
graph_df$continent_name_fac =  sortLvls.fnc(graph_df$continent_name_fac, 
                                            c(4, 3, 2, 1))

graph_df = graph_df %>% filter(year <= 2050)
graph_df$def_skm = graph_df$def_skm / 10
graph_df$clean_name_fac = factor(graph_df$clean_name)
graph_df$clean_name_fac_re = sortLvls.fnc(graph_df$clean_name_fac, 
                                          c(7,8, 5,6, 3,4, 1,2))

graph_df$year_1 = graph_df$year - 7.5
graph_df$year_2 = graph_df$year - 2.5
graph_df$def_skm_2 = graph_df$def_skm
graph_df$def_skm_2 = ((graph_df$def_skm_2 * 100) / 1000000) 
graph_df$def_skm = ((graph_df$def_skm * 100) / 1000000)

p <- ggplot() +
  geom_segment(data=graph_df, mapping=aes(x=year_1, y=def_skm, 
                                   xend=year_2, yend=def_skm_2, 
                                   color = continent_name_fac,
                                   linetype = clean_name_fac_re),
               size = 1.4) + 
  scale_linetype_manual(values=c("solid", "twodash", "solid", "twodash", "solid", "twodash", "solid", "twodash"),
                        guide = FALSE) +
  xlab('Year') +
  ylab("Deforestation (mha/yr)") +
  scale_x_continuous(breaks = seq(2000, 2050, by = 10),
                     limits = c(2000,2050)) +
  scale_y_continuous(limits = c(0, max(graph_df$def_skm)),
                     breaks = seq(0, max(graph_df$def_skm), by = 4),
                     labels = comma) +
  scale_colour_manual(values = brewer.pal(11, "BrBG")[c(1,2,3,4)], name = "") +
  #scale_colour_manual(values = brewer.pal(8, "Paired")[c(2,4,6,8)], name = "") +
  #ggtitle("Historical and projected tropical forest loss") +
  theme(
    plot.title=element_text(face="bold", size=13, colour="black")
    ,axis.text=element_text(size=11, colour = "black")
    ,axis.title=element_text(size=12,face="bold", colour = "black")
    ,legend.text = element_text(colour="black", size = 12)
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank())



png(paste0(graphics_folder_loc, "hist_ref_deforstation.png"), 
    width=11000,height=7600,res=1250, units="px")
print(p)
dev.off()


##### 1b. Create predicted business as usual deforestation
################################################################

graphics_folder_loc = "D:/Research/global_refor_mac/graphics/"

require(ggplot2)
require(scales)
require(RColorBrewer)

# match names 
match_df = data.frame("continent" = c("Africa",
                                      "Asia",
                                      "Latin America",
                                      "Pan-tropical"),
                      "name" = c("Africa", "Asia",
                                 "LatinAmerica", "global"),
                      stringsAsFactors = FALSE)

pred_emiss = standard_output[[3]] %>% filter(dollar == 0)
global_emiss = pred_emiss %>% group_by(year) %>% summarize(def_emiss = sum(def_emiss,
                                                                     na.rm = TRUE)) %>%
  mutate(continent = "global",
         dollar = 0)

hist_emiss = hist_data %>% group_by(cont_var) %>% summarize(def_emiss = sum(def_emiss,
                                                                        na.rm = TRUE))  %>%
  mutate(year = 2010,
         dollar = 0)
global_hist = hist_emiss %>% summarize(def_emiss = sum(def_emiss),
                                     year = mean(year),
                                     dollar = mean(dollar)) %>%
  mutate(cont_var = "global")

# create a historical data.frame 
hist_emiss = bind_rows(hist_emiss, global_hist)
pred_emiss = bind_rows(pred_emiss, global_emiss)
pred_emiss$year = as.numeric(pred_emiss$year)
colnames(hist_emiss)[1] = "continent"

# match the names 
hist_emiss$clean_name = paste0(match_df$continent[match(hist_emiss$continent, match_df$name)],
                             " - historical")
pred_emiss$clean_name = paste0(match_df$continent[match(pred_emiss$continent, match_df$name)],
                             " - projected")

# graph data.frame
graph_df = bind_rows(hist_emiss, pred_emiss)
graph_df = graph_df %>% filter(year <= 2050)
graph_df$def_emiss = graph_df$def_emiss / 10
graph_df$def_emiss = graph_df$def_emiss / 1000000000
graph_df$clean_name_fac = factor(graph_df$clean_name)
graph_df$clean_name_fac_re = sortLvls.fnc(graph_df$clean_name_fac, 
                                          c(7,8, 5,6, 3,4, 1,2))

graph_df$year_1 = graph_df$year - 7.5
graph_df$year_2 = graph_df$year - 2.5
graph_df$def_emiss_2 = graph_df$def_emiss

graph_df$continent_name = match_df$continent[match(graph_df$continent, match_df$name)]
graph_df$continent_name_fac = factor(graph_df$continent_name)
graph_df$continent_name_fac =  sortLvls.fnc(graph_df$continent_name_fac, 
                                            c(4, 3, 2, 1))

p <- ggplot() +
  geom_segment(data=graph_df, mapping=aes(x=year_1, y=def_emiss, 
                                          xend=year_2, yend=def_emiss_2, 
                                          color = continent_name_fac,
                                          linetype = clean_name_fac_re),
               size = 1.4) + 
  scale_linetype_manual(values=c("solid", "twodash", "solid", "twodash", "solid", "twodash", "solid", "twodash"),
                        guide = FALSE) +
  xlab('Year') +
  #ylab("Emissions from forest loss (GtCO2/yr)") +
  ylab(expression(paste("Emissions from deforestation (Gt", CO[2], "/yr)", sep=""))) +
  scale_x_continuous(breaks = seq(2000, 2050, by = 10),
                     limits = c(2000,2050)) +
  scale_y_continuous(limits = c(0, max(graph_df$def_emiss)),
                     breaks = seq(0, max(graph_df$def_emiss), by = 1),
                     labels = comma) +
  scale_colour_manual(values = brewer.pal(11, "BrBG")[c(1,2,3,4)], name = "") +
  #ggtitle("Historical and projected tropical forest loss") +
  theme(
    plot.title=element_text(face="bold", size=13, colour="black")
    ,axis.text=element_text(size=11, colour = "black")
    ,axis.title=element_text(size=12,face="bold", colour = "black")
    ,legend.text = element_text(colour="black", size = 12)
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank())



png(paste0(graphics_folder_loc, "hist_ref_deforstation_emissions.png"), 
    width=11000,height=7600,res=1250, units="px")
print(p)
dev.off()


##### 2 Create predicted business as usual reforestation
################################################################

graphics_folder_loc = "D:/Research/global_refor_mac/graphics/"

require(ggplot2)
require(scales)
require(RColorBrewer)

# match names 
match_df = data.frame("continent" = c("Africa",
                                      "Asia",
                                      "Latin America",
                                      "Pan-tropical"),
                      "name" = c("Africa", "Asia",
                                 "LatinAmerica", "global"),
                      stringsAsFactors = FALSE)

pred_ref = standard_output[[5]] %>% filter(dollar == 0)
global_ref = pred_ref %>% group_by(year) %>% summarize(gain_skm = sum(gain_skm,
                                                                      na.rm = TRUE)) %>%
  mutate(continent = "global",
         dollar = 0)

hist_ref = hist_data %>% group_by(cont_var) %>% summarize(gain_skm = sum(gain_skm,
                                                                         na.rm = TRUE))  %>%
  mutate(year = 2010,
         dollar = 0)
global_hist = hist_ref %>% summarize(gain_skm = sum(gain_skm),
                                     year = mean(year),
                                     dollar = mean(dollar)) %>%
  mutate(cont_var = "global")

# create a historical data.frame 
hist_ref = bind_rows(hist_ref, global_hist)
pred_ref = bind_rows(pred_ref, global_ref)
pred_ref$year = as.numeric(pred_ref$year)
colnames(hist_ref)[1] = "continent"

# match the names 
hist_ref$clean_name = paste0(match_df$continent[match(hist_ref$continent, match_df$name)],
                             " - historical")
pred_ref$clean_name = paste0(match_df$continent[match(pred_ref$continent, match_df$name)],
                             " - projected")

# graph data.frame
graph_df = bind_rows(hist_ref, pred_ref)
graph_df = graph_df %>% filter(year <= 2050)
graph_df.identifier = graph_df$clean_name

graph_df$gain_skm = graph_df$gain_skm / 10

graph_df$clean_name_fac = factor(graph_df$clean_name)
graph_df$clean_name_fac_re = sortLvls.fnc(graph_df$clean_name_fac, 
                                          c(7,8, 5,6, 3,4, 1,2))

graph_df$year_1 = graph_df$year - 7.5
graph_df$year_2 = graph_df$year - 2.5
graph_df$gain_skm_2 = graph_df$gain_skm

graph_df$continent_name = match_df$continent[match(graph_df$continent, match_df$name)]
graph_df$continent_name_fac = factor(graph_df$continent_name)
graph_df$continent_name_fac =  sortLvls.fnc(graph_df$continent_name_fac, 
                                            c(4, 3, 2, 1))


graph_df$gain_skm_2 = ((graph_df$gain_skm_2 * 100) / 1000000) 
graph_df$gain_skm = ((graph_df$gain_skm * 100) / 1000000)


p <- ggplot() +
  geom_segment(data=graph_df, mapping=aes(x=year_1, y=gain_skm, 
                                          xend=year_2, yend=gain_skm_2, 
                                          color = continent_name_fac,
                                          linetype = clean_name_fac_re),
               size = 1.4) + 
  xlab('Year') +
  ylab("Reforestation (mha/yr)") +
  scale_x_continuous(breaks = seq(2000, 2050, by = 10),
                     limits = c(2000,2050)) +
  scale_y_continuous(limits = c(0, max(graph_df$gain_skm)),
                     breaks = seq(0, max(graph_df$gain_skm), by = 2),
                     labels = comma) +
  scale_colour_manual(values = brewer.pal(9, "Greens")[c(9,7,5,3)], name = "") +
  #scale_colour_manual(values = brewer.pal(8, "Paired")[c(2,4,6,8)], name = "") +
  scale_linetype_manual(values=c("solid", "dashed", "solid", "twodash", "solid", "twodash", "solid", "twodash"), 
                        guide = F) +
  #ggtitle("Historical and projected tropical forest loss") +
  theme(
    plot.title=element_text(face="bold", size=13, colour="black")
    ,axis.text=element_text(size=11, colour = "black")
    ,axis.title=element_text(size=12,face="bold", colour = "black")
    ,legend.text = element_text(colour="black", size = 12)
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank())


png(paste0(graphics_folder_loc, "hist_ref_reforestation.png"), 
    width=11000,height=7600,res=1250, units="px")
print(p)
dev.off()


##### 2b Create predicted business as usual reforestation
################################################################

graphics_folder_loc = "D:/Research/global_refor_mac/graphics/"

require(ggplot2)
require(scales)
require(RColorBrewer)

# match names 
match_df = data.frame("continent" = c("Africa",
                                      "Asia",
                                      "Latin America",
                                      "Pan-tropical"),
                      "name" = c("Africa", "Asia",
                                 "LatinAmerica", "global"),
                      stringsAsFactors = FALSE)

pred_ref = standard_output[[6]] %>% filter(dollar == 0)
global_ref = pred_ref %>% group_by(year) %>% summarize(gain_emiss = sum(gain_emiss,
                                                                      na.rm = TRUE)) %>%
  mutate(continent = "global",
         dollar = 0)

hist_ref = hist_data %>% group_by(cont_var) %>% summarize(gain_emiss = sum(gain_emiss,
                                                                         na.rm = TRUE))  %>%
  mutate(year = 2010,
         dollar = 0)
global_hist = hist_ref %>% summarize(gain_emiss = sum(gain_emiss),
                                     year = mean(year),
                                     dollar = mean(dollar)) %>%
  mutate(cont_var = "global")

# create a historical data.frame 
hist_ref = bind_rows(hist_ref, global_hist)
pred_ref = bind_rows(pred_ref, global_ref)
pred_ref$year = as.numeric(pred_ref$year)
colnames(hist_ref)[1] = "continent"

# match the names 
hist_ref$clean_name = paste0(match_df$continent[match(hist_ref$continent, match_df$name)],
                             " - historical")
pred_ref$clean_name = paste0(match_df$continent[match(pred_ref$continent, match_df$name)],
                             " - projected")

# graph data.frame
graph_df = bind_rows(hist_ref, pred_ref)
graph_df = graph_df %>% filter(year <= 2050)
graph_df.identifier = graph_df$clean_name

graph_df$gain_emiss = graph_df$gain_emiss / 10
graph_df$gain_emiss = graph_df$gain_emiss / 1000000000

graph_df$clean_name_fac = factor(graph_df$clean_name)
graph_df$clean_name_fac_re = sortLvls.fnc(graph_df$clean_name_fac, 
                                          c(7,8, 5,6, 3,4, 1,2))

graph_df$year_1 = graph_df$year - 7.5
graph_df$year_2 = graph_df$year - 2.5
graph_df$gain_emiss_2 = graph_df$gain_emiss

graph_df$continent_name = match_df$continent[match(graph_df$continent, match_df$name)]
graph_df$continent_name_fac = factor(graph_df$continent_name)
graph_df$continent_name_fac =  sortLvls.fnc(graph_df$continent_name_fac, 
                                            c(4, 3, 2, 1))

p <- ggplot() +
  geom_segment(data=graph_df, mapping=aes(x=year_1, y=gain_emiss, 
                                          xend=year_2, yend=gain_emiss_2, 
                                          color = continent_name_fac,
                                          linetype = clean_name_fac_re),
               size = 1.4) + 
  xlab('Year') +
  #ylab("Removals from forest gain (GtCO2/yr)") +
  ylab(expression(paste("Removals from reforestation (Gt", CO[2], "/yr)", sep=""))) +
  scale_x_continuous(breaks = seq(2000, 2050, by = 10),
                     limits = c(2000,2050)) +
  scale_y_continuous(limits = c(0, max(graph_df$gain_emiss)),
                     breaks = seq(0, max(graph_df$gain_emiss), by = 0.5),
                     labels = comma) +
  scale_colour_manual(values = brewer.pal(9, "Greens")[c(9,7,5,3)], name = "") +
  scale_linetype_manual(values=c("solid", "dashed", "solid", "twodash", "solid", "twodash", "solid", "twodash"),
                        guide = FALSE) +
  #ggtitle("Historical and projected tropical forest loss") +
  theme(
    plot.title=element_text(face="bold", size=13, colour="black")
    ,axis.text=element_text(size=11, colour = "black")
    ,axis.title=element_text(size=12,face="bold", colour = "black")
    ,legend.text = element_text(colour="black", size = 12)
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank())


png(paste0(graphics_folder_loc, "hist_ref_reforestation_emissions.png"), 
    width=11000,height=7600,res=1250, units="px")
print(p)
dev.off()


##### 3 Calculate projected deforestation per year under BAU w carbon prices
################################################################

graphics_folder_loc = "D:/Research/global_refor_mac/graphics/"

require(ggplot2)
require(scales)
require(RColorBrewer)

# match names 
match_df = data.frame("dollar" = c(0, 10, 20, 50, 100),
                      "name" = c("BAU", "$10",
                                 "$20", "$50", "$100"),
                      stringsAsFactors = FALSE)

pred_def_skm = standard_output[[2]] %>% filter(dollar %in% c(0, 10, 20, 50, 100)) %>% 
  group_by(year,dollar) %>%
  summarize(def_skm = sum(def_skm, na.rm = TRUE)) 

hist_def_skm = hist_data %>% summarize(def_skm = sum(def_skm,
                                                         na.rm = TRUE))  %>%
  summarize(def_skm = sum(def_skm, na.rm = TRUE)) %>%
  mutate(year = 2010,
         name_cat = "Historical") 

first_def_skm = pred_def_skm %>% filter(year == 2020) %>%
          mutate(name_cat = "BAU")
pred_def_skm = pred_def_skm %>% filter(year > 2020)
  

# bind the rows
pred_def_skm$year = as.numeric(pred_def_skm$year)
hist_def_skm$year = as.numeric(hist_def_skm$year)
first_def_skm$year = as.numeric(first_def_skm$year)
pred_def_skm$name_cat = match_df$name[match(pred_def_skm$dollar, match_df$dollar)]

graph_skm = bind_rows(first_def_skm, hist_def_skm)
graph_skm = bind_rows(graph_skm, pred_def_skm)
graph_skm$def_skm = graph_skm$def_skm / 10
graph_skm = graph_skm %>% filter(year <= 2050)

require(ggplot2)
require(scales)
require(RColorBrewer)

graph_skm$name_cat_fac = factor(graph_skm$name_cat)
graph_skm$name_cat_fac_re = sortLvls.fnc(graph_skm$name_cat_fac, 
                                          c(6, 5, 1, 3, 4, 2))

graph_skm$year_1 = graph_skm$year - 7.5
graph_skm$year_2 = graph_skm$year - 2.5
graph_skm$def_skm_2 = graph_skm$def_skm

graph_skm$def_skm_2 = ((graph_skm$def_skm_2 * 100) / 1000000) 
graph_skm$def_skm = ((graph_skm$def_skm * 100) / 1000000)

g1 <- ggplot() +
  geom_segment(data=graph_skm, mapping=aes(x=year_1, y=def_skm, 
                                          xend=year_2, yend=def_skm_2, 
                                          color = name_cat_fac_re,
                                          linetype = name_cat_fac_re),
               size = 1.4) + 
  xlab('Year') +
  ylab('Deforestation (mha/yr)') +
  scale_x_continuous(breaks = seq(2000, 2050, by = 10),
                     limits = c(2000,2050)) +  
  scale_y_continuous(limits = c(0, max(graph_skm$def_skm)),
                     breaks = seq(0, max(graph_skm$def_skm), by = 4),
                     label = comma) +
  scale_linetype_manual(values=c("solid", "twodash", "twodash", "twodash", "twodash", 
                                 "twodash", "twodash")) +
  scale_colour_manual(values = c("#281600", "#5b3507", brewer.pal(11, "BrBG")[c(2,3,4, 5)]), name = "") +
  # scale_colour_manual(values = brewer.pal(6, "Spectral")) +
  #ggtitle("Historical and projected emissions from tropical forest loss, by carbon price ($/tCO2)") +
  theme(
    plot.title=element_text(face="bold", size=12, colour="black")
    ,axis.text.y=element_text(size=10, colour = "black")
    ,axis.text.x=element_text(size=9, colour = "black")
    ,axis.title=element_text(size=11,face="bold", colour = "black")
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank()
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,legend.text = element_text(colour="black", size = 11, face = "bold")
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
  ) 

png(paste0(graphics_folder_loc, "projected_deforestation_skm.png"), 
    width=11000,height=7600,res=1250, units="px")
print(g1)
dev.off()


##### 4 Calculate projected deforestation per year under BAU w carbon prices
################################################################

graphics_folder_loc = "D:/Research/global_refor_mac/graphics/"

require(ggplot2)
require(scales)
require(RColorBrewer)

# match names 
match_df = data.frame("dollar" = c(0, 10, 20, 50, 100),
                      "name" = c("BAU", "$10",
                                 "$20", "$50", "$100"),
                      stringsAsFactors = FALSE)

pred_def_emiss = standard_output[[3]] %>% filter(dollar %in% c(0, 10, 20, 50, 100)) %>% 
  group_by(year,dollar) %>%
  summarize(def_emiss = sum(def_emiss, na.rm = TRUE)) 

hist_def_emiss = hist_data %>% summarize(def_emiss = sum(def_emiss,
                                                         na.rm = TRUE))  %>%
  summarize(def_emiss = sum(def_emiss, na.rm = TRUE)) %>%
  mutate(year = 2010,
         name_cat = "Historical") 

first_def_emiss = pred_def_emiss %>% filter(year == 2020) %>%
  mutate(name_cat = "BAU")
pred_def_emiss = pred_def_emiss %>% filter(year > 2020)


# bind the rows
pred_def_emiss$year = as.numeric(pred_def_emiss$year)
hist_def_emiss$year = as.numeric(hist_def_emiss$year)
first_def_emiss$year = as.numeric(first_def_emiss$year)
pred_def_emiss$name_cat = match_df$name[match(pred_def_emiss$dollar, match_df$dollar)]

graph_emiss = bind_rows(first_def_emiss, hist_def_emiss)
graph_emiss = bind_rows(graph_emiss, pred_def_emiss)
graph_emiss$def_emiss = graph_emiss$def_emiss / 10
graph_emiss$def_emiss = graph_emiss$def_emiss / 1000000000
graph_emiss = graph_emiss %>% filter(year <= 2050)

require(ggplot2)
require(scales)
require(RColorBrewer)


graph_emiss$year_1 = graph_emiss$year - 7.5
graph_emiss$year_2 = graph_emiss$year - 2.5
graph_emiss$def_emiss_2 = graph_emiss$def_emiss

graph_emiss$name_cat_fac = factor(graph_emiss$name_cat)
graph_emiss$name_cat_fac_re = sortLvls.fnc(graph_emiss$name_cat_fac, 
                                         c(6, 5, 1, 3, 4, 2))


g1 <- ggplot() +
  geom_segment(data=graph_emiss, mapping=aes(x=year_1, y=def_emiss, 
                                           xend=year_2, yend=def_emiss_2, 
                                           color = name_cat_fac_re,
                                           linetype = name_cat_fac_re),
               size = 1.4) + 
  xlab('Year') +
  ylab(expression(paste("Emissions from deforestation (Gt", CO[2], "/yr)", sep=""))) +
  scale_x_continuous(breaks = seq(2000, 2050, by = 10),
                     limits = c(2000,2050)) +  
  scale_y_continuous(limits = c(0, max(graph_emiss$def_emiss)),
                     breaks = seq(0, max(graph_emiss$def_emiss), by = 1),
                     label = comma) +
  scale_linetype_manual(values=c("solid", "twodash", "twodash", "twodash", "twodash", 
                                 "twodash"), name = "") +
  scale_colour_manual(values = c("#281600", "#5b3507", brewer.pal(11, "BrBG")[c(2,3,4, 5)]), name = "") +
  #ggtitle("Historical and projected emissions from tropical forest loss, by carbon price ($/tCO2)") +
  theme(
    plot.title=element_text(face="bold", size=12, colour="black")
    ,axis.text.y=element_text(size=10, colour = "black")
    ,axis.text.x=element_text(size=9, colour = "black")
    ,axis.title=element_text(size=11,face="bold", colour = "black")
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank()
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,legend.text = element_text(colour="black", size = 11, face = "bold")
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
  ) 


png(paste0(graphics_folder_loc, "projected_deforestation_emissions.png"), 
    width=11000,height=7600,res=1250, units="px")
print(g1)
dev.off()




##### 5 Calculate projected reforestation per year under BAU w carbon prices
################################################################

graphics_folder_loc = "D:/Research/global_refor_mac/graphics/"

require(ggplot2)
require(scales)
require(RColorBrewer)

# match names 
match_df = data.frame("dollar" = c(0, 10, 20, 50, 100),
                      "name" = c("BAU", "$10",
                                 "$20", "$50", "$100"),
                      stringsAsFactors = FALSE)

pred_gain_skm = standard_output[[5]] %>% filter(dollar %in% c(0, 10, 20, 50, 100)) %>% 
  group_by(year,dollar) %>%
  summarize(gain_skm = sum(gain_skm, na.rm = TRUE)) 

hist_gain_skm = hist_data %>% summarize(gain_skm = sum(gain_skm,
                                                     na.rm = TRUE))  %>%
  summarize(gain_skm = sum(gain_skm, na.rm = TRUE)) %>%
  mutate(year = 2010,
         name_cat = "Historical") 

first_gain_skm = pred_gain_skm %>% filter(year == 2020) %>%
  mutate(name_cat = "BAU")
pred_gain_skm = pred_gain_skm %>% filter(year > 2020)

# bind the rows
pred_gain_skm$year = as.numeric(pred_gain_skm$year)
hist_gain_skm$year = as.numeric(hist_gain_skm$year)
first_gain_skm$year = as.numeric(first_gain_skm$year)
pred_gain_skm$name_cat = match_df$name[match(pred_gain_skm$dollar, match_df$dollar)]

graph_skm = bind_rows(first_gain_skm, hist_gain_skm)
graph_skm = bind_rows(pred_gain_skm, graph_skm)
graph_skm$gain_skm = graph_skm$gain_skm / 10
graph_skm = graph_skm %>% filter(year <= 2050)

require(ggplot2)
require(scales)
require(RColorBrewer)

graph_skm$name_cat_fac = factor(graph_skm$name_cat)
graph_skm$name_cat_fac_re = sortLvls.fnc(graph_skm$name_cat_fac, 
                                           c(6, 5, 1, 3, 4, 2))


graph_skm$year_1 = graph_skm$year - 7.5
graph_skm$year_2 = graph_skm$year - 2.5
graph_skm$gain_skm_2 = graph_skm$gain_skm

graph_skm$gain_skm_2 = ((graph_skm$gain_skm_2 * 100) / 1000000) 
graph_skm$gain_skm = ((graph_skm$gain_skm * 100) / 1000000)


g1 <-ggplot() +
  geom_segment(data=graph_skm, mapping=aes(x=year_1, y=gain_skm, 
                                             xend=year_2, yend=gain_skm_2, 
                                             color = name_cat_fac_re,
                                             linetype = name_cat_fac_re),
               size = 1.4) + 
  xlab('Year') +
  ylab('Reforestation (mha/yr)') +
  scale_x_continuous(breaks = seq(2000, 2050, by = 10),
                     limits = c(2000,2050)) +  
  scale_y_continuous(limits = c(0, max(graph_skm$gain_skm)),
                     breaks = seq(0, max(graph_skm$gain_skm), by = 3),
                     label = comma) +
  scale_linetype_manual(values=c("solid", "twodash", "twodash", "twodash", "twodash", 
                                 "twodash"), name = "") +
  scale_colour_manual(values = c("#001c0b", brewer.pal(9, "Greens")[c(9,7,5,3, 2)]), name = "") +
  #ggtitle("Historical and projected emissions from tropical forest loss, by carbon price ($/tCO2)") +
  theme(
    plot.title=element_text(face="bold", size=12, colour="black")
    ,axis.text.y=element_text(size=10, colour = "black")
    ,axis.text.x=element_text(size=9, colour = "black")
    ,axis.title=element_text(size=11,face="bold", colour = "black")
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank()
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,legend.text = element_text(colour="black", size = 11, face = "bold")
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
  ) 

png(paste0(graphics_folder_loc, "projected_reforestation_skm.png"), 
    width=11000,height=7600,res=1250, units="px")
print(g1)
dev.off()




##### 6 Calculate projected reforestation emissions per year under BAU w carbon prices
################################################################

graphics_folder_loc = "D:/Research/global_refor_mac/graphics/"

require(ggplot2)
require(scales)
require(RColorBrewer)

# match names 
match_df = data.frame("dollar" = c(0, 10, 20, 50, 100),
                      "name" = c("BAU", "$10",
                                 "$20", "$50", "$100"),
                      stringsAsFactors = FALSE)

pred_def_emiss = standard_output[[6]] %>% filter(dollar %in% c(0, 10, 20, 50, 100)) %>% 
  group_by(year,dollar) %>%
  summarize(gain_emiss = sum(gain_emiss, na.rm = TRUE)) 

hist_def_emiss = hist_data %>% summarize(gain_emiss = sum(gain_emiss,
                                                         na.rm = TRUE))  %>%
  summarize(gain_emiss = sum(gain_emiss, na.rm = TRUE)) %>%
  mutate(year = 2010,
         name_cat = "Historical") 

first_gain_emiss = pred_def_emiss %>% filter(year == 2020) %>%
  mutate(name_cat = "BAU")
pred_def_emiss = pred_def_emiss %>% filter(year > 2020)

# bind the rows
pred_def_emiss$year = as.numeric(pred_def_emiss$year)
hist_def_emiss$year = as.numeric(hist_def_emiss$year)
first_gain_emiss$year = as.numeric(first_gain_emiss$year)
pred_def_emiss$name_cat = match_df$name[match(pred_def_emiss$dollar, match_df$dollar)]

graph_emiss = bind_rows(first_gain_emiss, hist_def_emiss)
graph_emiss = bind_rows(graph_emiss, pred_def_emiss)
graph_emiss$gain_emiss = graph_emiss$gain_emiss / 10
graph_emiss$gain_emiss = graph_emiss$gain_emiss / 1000000000
graph_emiss = graph_emiss %>% filter(year <= 2050)

require(ggplot2)
require(scales)
require(RColorBrewer)

graph_emiss$name_cat_fac = factor(graph_emiss$name_cat)
graph_emiss$name_cat_fac_re = sortLvls.fnc(graph_emiss$name_cat_fac, 
                                         c(6, 5, 1, 3, 4, 2))


graph_emiss$year_1 = graph_emiss$year - 7.5
graph_emiss$year_2 = graph_emiss$year - 2.5
graph_emiss$gain_emiss_2 = graph_emiss$gain_emiss


g1 <- ggplot() +
  geom_segment(data=graph_emiss, mapping=aes(x=year_1, y=gain_emiss, 
                                           xend=year_2, yend=gain_emiss_2, 
                                           color = name_cat_fac_re,
                                           linetype = name_cat_fac_re),
               size = 1.4) +  
  xlab('Year') +
  ylab(expression(paste("Removals from reforestation (Gt", CO[2], "/yr)", sep=""))) +
  scale_x_continuous(breaks = seq(2000, 2050, by = 10),
                     limits = c(2000,2050)) +  
  scale_y_continuous(limits = c(0, max(graph_emiss$gain_emiss)),
                     breaks = seq(0, max(graph_emiss$gain_emiss), by = 0.5),
                     label = comma) +
  scale_linetype_manual(values=c("solid", "twodash", "twodash", "twodash", "twodash", 
                                 "twodash"), name = "") +
  scale_colour_manual(values = c("#001c0b", brewer.pal(9, "Greens")[c(9,7,5,3, 2)]), name = "") +
  #ggtitle("Historical and projected emissions from tropical forest loss, by carbon price ($/tCO2)") +
  theme(
    plot.title=element_text(face="bold", size=12, colour="black")
    ,axis.text.y=element_text(size=10, colour = "black")
    ,axis.text.x=element_text(size=9, colour = "black")
    ,axis.title=element_text(size=11,face="bold", colour = "black")
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank()
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,legend.text = element_text(colour="black", size = 11, face = "bold")
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
  ) 


png(paste0(graphics_folder_loc, "projected_reforestation_emissions.png"), 
    width=11000,height=7600,res=1250, units="px")
print(g1)
dev.off()


######################################################################
defor_for_skm = standard_output[[1]]
defor_for_skm$fc_skm = defor_for_skm$fc_skm / 10
defor_for_skm = defor_for_skm %>% spread(year, fc_skm)

def_skm = standard_output[[2]]
def_skm$def_skm = def_skm$def_skm / 10
def_skm = def_skm %>% spread(year, def_skm)

gain_skm = standard_output[[5]]
gain_skm$gain_skm = gain_skm$gain_skm / 10
gain_skm = gain_skm %>% spread(year, gain_skm)

refor_for_skm = standard_output[[4]]
refor_for_skm$fc_skm = refor_for_skm$fc_skm / 10
refor_for_skm = refor_for_skm %>% spread(year, fc_skm)

def_emiss = standard_output[[3]]
def_emiss$def_emiss = def_emiss$def_emiss / 10
def_emiss$def_emiss = def_emiss$def_emiss / 1000000000
def_emiss = def_emiss %>% spread(year, def_emiss)

gain_emiss = standard_output[[6]]
gain_emiss$gain_emiss = gain_emiss$gain_emiss / 10
gain_emiss$gain_emiss = gain_emiss$gain_emiss / 1000000000
gain_emiss = gain_emiss %>% spread(year, gain_emiss)

write.csv(defor_for_skm, paste0("D:/Research/global_refor_mac/graphics/global_forest_cover_skm_only_defor_model.csv"))
write.csv(def_skm, paste0("D:/Research/global_refor_mac/graphics/global_deforestation_skm.csv"))
write.csv(gain_skm, paste0("D:/Research/global_refor_mac/graphics/global_reforestation_skm.csv"))
write.csv(refor_for_skm, paste0("D:/Research/global_refor_mac/graphics/global_forest_cover_skm_only_refor_model.csv"))
write.csv(def_emiss, paste0("D:/Research/global_refor_mac/graphics/global_deforestation_emissions_gigatonnes.csv"))
write.csv(gain_emiss, paste0("D:/Research/global_refor_mac/graphics/global_reforestation_removals_gigatonnes.csv"))


##########################################################################################
### def_emiss = standard_output[[3]]
## MAC deforestation co2 emissions
def_emiss = standard_output[[3]] %>%
                group_by(continent, year) %>%
                arrange(continent, year, dollar) %>%
                mutate(mac_emiss = abs(def_emiss - first(def_emiss)) / 10)
def_emiss = def_emiss %>% filter(year <= 2050)
def_emiss$mac_emiss = def_emiss$mac_emiss / 1000000000

glob_def_emiss = def_emiss %>% 
                  group_by(year, dollar) %>%
                  summarize(mac_emiss = sum(mac_emiss),
                            def_emiss = sum(def_emiss)) %>%
                  mutate(continent = "global")
glob_def_emiss = glob_def_emiss %>% filter(year <= 2050) %>% filter(year > 2020)

# bind rows
def_emiss = bind_rows(def_emiss, glob_def_emiss)

g1 <- ggplot(glob_def_emiss, mapping=aes(x = mac_emiss, y = dollar, 
                                         color = year, group = year)) +
  geom_line(size = 1.4) +
  xlab('Reduced emissions (GtCO2/yr)') +
  geom_vline(xintercept=seq(0, ceiling(max(glob_def_emiss$mac_emiss)), by = 1), linetype="dotted") + 
  ylab('Carbon price ($/t CO2)') +
  scale_y_continuous(breaks = seq(0, 100, by = 10),
                     limits = c(0,100)) +  
  scale_x_continuous(limits = c(0, ceiling(max(glob_def_emiss$mac_emiss))),
                     breaks = seq(0, ceiling(max(glob_def_emiss$mac_emiss)), by = 1),
                     label = comma) +
  scale_colour_manual(values = brewer.pal(3, "Spectral"), name = "", label = c("2021 - 2030",
                                                                               "2031 - 2040",
                                                                               "2041 - 2050")) +
  #ggtitle("Historical and projected emissions from tropical forest loss, by carbon price ($/tCO2)") +
  theme(
    plot.title=element_text(face="bold", size=12, colour="black")
    ,axis.text.y=element_text(size=10, colour = "black")
    ,axis.text.x=element_text(size=9, colour = "black")
    ,axis.title=element_text(size=11,face="bold", colour = "black")
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank()
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,legend.text = element_text(colour="black", size = 11, face = "bold")
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
  ) 


png(paste0(graphics_folder_loc, "mac_deforestation_co2.png"), 
    width=11000,height=7600,res=1250, units="px")
print(g1)
dev.off()




## MAC reforestation co2 emissions
ref_emiss = standard_output[[6]] %>%
  group_by(continent, year) %>%
  arrange(continent, year, dollar) %>%
  mutate(mac_emiss = abs(gain_emiss - first(gain_emiss)) / 10)
ref_emiss = ref_emiss %>% filter(year <= 2050)
ref_emiss$mac_emiss = ref_emiss$mac_emiss / 1000000000

glob_ref_emiss = ref_emiss %>% 
  group_by(year, dollar) %>%
  summarize(mac_emiss = sum(mac_emiss),
            gain_emiss = sum(gain_emiss)) %>%
  mutate(continent = "global")
glob_ref_emiss = glob_ref_emiss %>% filter(year <= 2050) %>% filter(year > 2020)
# glob_ref_emiss$mac_emiss = glob_ref_emiss$mac_emiss * 1000

# bind rows
ref_emiss = bind_rows(ref_emiss, glob_ref_emiss)

g1 <- ggplot(glob_ref_emiss, mapping=aes(x = mac_emiss, y = dollar, 
                                         color = year, group = year)) +
  geom_line(size = 1.4) +
  xlab('Enhanced removals (GtCO2/yr)') +
  geom_vline(xintercept=seq(0, ceiling(max(glob_def_emiss$mac_emiss)), by = 1), 
             linetype="dotted") + 
  ylab('Carbon price ($/t CO2)') +
  scale_y_continuous(breaks = seq(0, 100, by = 10),
                     limits = c(0,100)) +  
  scale_x_continuous(limits = c(0, ceiling(max(glob_def_emiss$mac_emiss))),
                     breaks = seq(0, ceiling(max(glob_def_emiss$mac_emiss)), by = 1),
                     label = comma) +
  scale_colour_manual(values = brewer.pal(3, "Spectral"), name = "", label = c( "2021 - 2030",
                                                                               "2031 - 2040",
                                                                               "2041 - 2050")) +
  #ggtitle("Historical and projected emissions from tropical forest loss, by carbon price ($/tCO2)") +
  theme(
    plot.title=element_text(face="bold", size=12, colour="black")
    ,axis.text.y=element_text(size=10, colour = "black")
    ,axis.text.x=element_text(size=9, colour = "black")
    ,axis.title=element_text(size=11,face="bold", colour = "black")
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank()
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,legend.text = element_text(colour="black", size = 11, face = "bold")
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
  ) 



png(paste0(graphics_folder_loc, "mac_reforestation_co2.png"), 
    width=11000,height=7600,res=1250, units="px")
print(g1)
dev.off()

# bind dataset
glob_def_emiss$category = "Avoided deforestation"
glob_ref_emiss$category = "Enhanced reforestation"
glob_df = bind_rows(glob_def_emiss, glob_ref_emiss)
glob_df$full_category = paste0(glob_df$category, ', ', as.numeric(glob_df$year) - 9, " - ",  as.numeric(glob_df$year))
glob_df$full_category_fac = factor(glob_df$full_category )
glob_df$full_category_fac = sortLvls.fnc(glob_df$full_category_fac, 
                                           c(4, 5, 6, 1, 2, 3))


g1 <- ggplot(glob_df, mapping=aes(x = mac_emiss, y = dollar, 
                                         color = full_category_fac, group = full_category_fac)) +
  geom_line(size = 1.43) +
  #xlab('Abatement (GtCO2/yr)') +
  #geom_vline(xintercept=seq(0, ceiling(max(glob_def_emiss$mac_emiss)), by = 1), 
  #           linetype="dotted") + 
  xlab(expression(paste("Abatement (Gt", CO[2], "/yr)", sep=""))) +
  ylab(expression(paste("Carbon price ($/t", CO[2], ")", sep=""))) +
  scale_y_continuous(breaks = seq(0, 100, by = 10),
                     limits = c(0,100),
                     expand = c(0,0)) +  
  scale_x_continuous(limits = c(0, ceiling(max(glob_def_emiss$mac_emiss))),
                     breaks = seq(0, ceiling(max(glob_def_emiss$mac_emiss)), by = 1),
                     label = comma,
                     expand = c(0,0)) +
  scale_colour_manual(values = c(brewer.pal(9, "Greens")[c(3,5,8)],
                                 c("#f6e8c3", "#bf812d", "#543005")), name = "",
                      label = c("Increased removals, \n 2021 - 2030",
                                "2031 - 2040",
                                "2041 - 2050",
                                "Reduced emissions, \n 2021 - 2030",
                                "2031 - 2040",
                                "2041 - 2050")) +
  #ggtitle("Historical and projected emissions from tropical forest loss, by carbon price ($/tCO2)") +
  theme(
    plot.title=element_text(face="bold", size=12, colour="black")
    ,axis.text.y=element_text(size=10, colour = "black")
    ,axis.text.x=element_text(size=9, colour = "black")
    ,axis.title=element_text(size=11,face="bold", colour = "black")
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.background = element_blank()
    ,legend.title = element_blank()
    ,axis.line.x = element_line(size = 1, colour = "black", linetype = 1)
    ,axis.line.y = element_line(size = 1, colour = "black", linetype = 1)
    ,legend.text = element_text(colour="black", size = 11, face = "bold")
    ,legend.key = element_rect(colour = "white", fill = "white")
    ,axis.ticks = element_line(colour = 'black')
    ,legend.key.size = unit(2, 'lines')
  ) 



png(paste0(graphics_folder_loc, "mac_deforestation_reforestation_co2.png"), 
    width=11000,height=7600,res=1250, units="px")
print(g1)
dev.off()

#############################################################
#### create mac curve for reforestation

ref_emiss = standard_output[[6]] %>%
  group_by(continent, year) %>%
  arrange(continent, year, dollar) %>%
  mutate(mac_emiss = abs(gain_emiss - first(gain_emiss)) / 10)
ref_emiss = ref_emiss %>% filter(year <= 2050)
ref_emiss$mac_emiss = ref_emiss$mac_emiss / 1000000000

glob_ref_emiss = ref_emiss %>% 
  group_by(year, dollar) %>%
  summarize(mac_emiss = sum(mac_emiss),
            gain_emiss = sum(gain_emiss)) %>%
  mutate(continent = "Global")
glob_ref_emiss = glob_ref_emiss %>% filter(year <= 2050) %>% filter(year > 2020)

# bind rows
ref_emiss = bind_rows(ref_emiss, glob_ref_emiss)

mac_ref_emiss = data.frame(ref_emiss %>% select(-gain_emiss) %>% spread(dollar, mac_emiss),
                           stringsAsFactors = FALSE)
colnames(mac_ref_emiss) = c("Continent", "Year", paste0("carbon_price_",
                                                        c(0, 5, seq(10, 100, 10))))

write.csv(mac_ref_emiss, paste0("D:/Research/global_refor_mac/graphics/mac_curve_by_continent_and_global_refor_gigatonnes.csv"))


## mac deforestation
def_emiss = standard_output[[3]] %>%
  group_by(continent, year) %>%
  arrange(continent, year, dollar) %>%
  mutate(mac_emiss = abs(def_emiss - first(def_emiss)) / 10)
def_emiss = def_emiss %>% filter(year <= 2050)
def_emiss$mac_emiss = def_emiss$mac_emiss / 1000000000

glob_def_emiss = def_emiss %>% 
  group_by(year, dollar) %>%
  summarize(mac_emiss = sum(mac_emiss),
            def_emiss = sum(def_emiss)) %>%
  mutate(continent = "Global")
glob_def_emiss = glob_def_emiss %>% filter(year <= 2050) %>% filter(year > 2020)

# bind rows
def_emiss = bind_rows(def_emiss, glob_def_emiss)

mac_def_emiss = data.frame(def_emiss %>% select(-def_emiss) %>% spread(dollar, mac_emiss),
                           stringsAsFactors = FALSE)
colnames(mac_def_emiss) = c("Continent", "Year", paste0("carbon_price_",
                                                        c(0, 5, seq(10, 100, 10))))

write.csv(mac_def_emiss, paste0("D:/Research/global_refor_mac/graphics/mac_curve_by_continent_and_global_defor_gigatonnes.csv"))



#############################################################################
standard_output_cntry <- calc_data_model_cntry(model_list = model_list,
                                               year_vec = seq(2020, 2050, 10),
                                               dollar_prices = c(0,5,seq(10,100,10)))

# calculate the reforestation amount
defor_skm_cntry = standard_output_cntry[[2]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 0) %>%
  group_by(country) %>%
  summarize(def_skm = sum(def_skm, na.rm = TRUE))

refor_skm_cntry = standard_output_cntry[[5]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 0) %>%
  group_by(country) %>%
  summarize(gain_skm = sum(gain_skm, na.rm = TRUE))

defor_skm_cntry_20dol = standard_output_cntry[[2]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 20) %>%
  group_by(country) %>%
  summarize(def_skm = sum(def_skm, na.rm = TRUE))

refor_skm_cntry_20dol = standard_output_cntry[[5]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 20) %>%
  group_by(country) %>%
  summarize(gain_skm = sum(gain_skm, na.rm = TRUE))

defor_skm_cntry_50dol = standard_output_cntry[[2]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 50) %>%
  group_by(country) %>%
  summarize(def_skm = sum(def_skm, na.rm = TRUE))

refor_skm_cntry_50dol = standard_output_cntry[[5]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 50) %>%
  group_by(country) %>%
  summarize(gain_skm = sum(gain_skm, na.rm = TRUE))


#### reforestation amounts
defor_emiss_cntry = standard_output_cntry[[3]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 0) %>%
  group_by(country) %>%
  summarize(def_emiss = sum(def_emiss, na.rm = TRUE))

refor_emiss_cntry = standard_output_cntry[[6]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 0) %>%
  group_by(country) %>%
  summarize(gain_emiss = sum(gain_emiss, na.rm = TRUE))

defor_emiss_cntry_20dol = standard_output_cntry[[3]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 20) %>%
  group_by(country) %>%
  summarize(def_emiss = sum(def_emiss, na.rm = TRUE))

refor_emiss_cntry_20dol = standard_output_cntry[[6]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 20) %>%
  group_by(country) %>%
  summarize(gain_emiss = sum(gain_emiss, na.rm = TRUE))

defor_emiss_cntry_50dol = standard_output_cntry[[3]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 50) %>%
  group_by(country) %>%
  summarize(def_emiss = sum(def_emiss, na.rm = TRUE))

refor_emiss_cntry_50dol = standard_output_cntry[[6]] %>% 
  filter(year %in% seq(2030, 2050, 10)) %>%
  filter(dollar == 50) %>%
  group_by(country) %>%
  summarize(gain_emiss = sum(gain_emiss, na.rm = TRUE))



### merge all of the deforestation data
colnames(defor_skm_cntry_20dol) = c("country", "def_skm_20dol")
colnames(defor_skm_cntry_50dol) = c("country", "def_skm_50dol")
colnames(defor_skm_cntry) = c("country", "def_skm_0dol")
defor_skm_cntry = left_join(defor_skm_cntry, defor_skm_cntry_20dol, by = "country")
defor_skm_cntry = left_join(defor_skm_cntry, defor_skm_cntry_50dol, by = "country")
defor_skm_cntry$red_def_skm_20dol = defor_skm_cntry$def_skm_0dol - defor_skm_cntry$def_skm_20dol
defor_skm_cntry$red_def_skm_50dol = defor_skm_cntry$def_skm_0dol - defor_skm_cntry$def_skm_50dol


### merge all of the deforestation emissions data
colnames(defor_emiss_cntry_20dol) = c("country", "def_emiss_20dol")
colnames(defor_emiss_cntry_50dol) = c("country", "def_emiss_50dol")
colnames(defor_emiss_cntry) = c("country", "def_emiss_0dol")
defor_emiss_cntry = left_join(defor_emiss_cntry, defor_emiss_cntry_20dol, by = "country")
defor_emiss_cntry = left_join(defor_emiss_cntry, defor_emiss_cntry_50dol, by = "country")
defor_emiss_cntry$avoid_emiss_def_20dol = defor_emiss_cntry$def_emiss_0dol - defor_emiss_cntry$def_emiss_20dol
defor_emiss_cntry$avoid_emiss_def_50dol = defor_emiss_cntry$def_emiss_0dol - defor_emiss_cntry$def_emiss_50dol


### merge all of the deforestation data
colnames(refor_skm_cntry_20dol) = c("country", "ref_skm_20dol")
colnames(refor_skm_cntry_50dol) = c("country", "ref_skm_50dol")
colnames(refor_skm_cntry) = c("country", "ref_skm_0dol")
refor_skm_cntry = left_join(refor_skm_cntry, refor_skm_cntry_20dol, by = "country")
refor_skm_cntry = left_join(refor_skm_cntry, refor_skm_cntry_50dol, by = "country")
refor_skm_cntry$enh_ref_skm_20dol = refor_skm_cntry$ref_skm_20dol - refor_skm_cntry$ref_skm_0dol
refor_skm_cntry$enh_ref_skm_50dol = refor_skm_cntry$ref_skm_50dol - refor_skm_cntry$ref_skm_0dol

### merge all of the reforestation data
colnames(refor_emiss_cntry_20dol) = c("country", "ref_emiss_20dol")
colnames(refor_emiss_cntry_50dol) = c("country", "ref_emiss_50dol")
colnames(refor_emiss_cntry) = c("country", "ref_emiss_0dol")
refor_emiss_cntry = left_join(refor_emiss_cntry, refor_emiss_cntry_20dol, by = "country")
refor_emiss_cntry = left_join(refor_emiss_cntry, refor_emiss_cntry_50dol, by = "country")
refor_emiss_cntry$enh_emiss_ref_20dol = refor_emiss_cntry$ref_emiss_20dol - refor_emiss_cntry$ref_emiss_0dol
refor_emiss_cntry$enh_emiss_ref_50dol = refor_emiss_cntry$ref_emiss_50dol - refor_emiss_cntry$ref_emiss_0dol


# merge dataset
full_out_df = left_join(defor_skm_cntry, defor_emiss_cntry, by = "country")
full_out_df = left_join(full_out_df, refor_skm_cntry, by = "country")
full_out_df = left_join(full_out_df, refor_emiss_cntry, by = "country")
full_out_df = full_out_df %>% arrange(desc(def_skm_0dol)) %>%
                          filter(def_skm_0dol > 0)



write.csv(full_out_df, 
          paste0("D:/Research/global_refor_mac/graphics/country_statistics_deforestation_reforestation.csv"))
